library(elevatr)
library(movecost)
library(sf)
library(rgdal)
library(leastcostpath)
library(sp)
library(raster)
library(tidyverse)
library(ggplot2)
library(dplyr)
library(terra)
library(ncdf4) 
library(exactextractr)  
library(rgeos)
library(PBSmapping)
library(geosphere)

setwd("/Users/Martha/Dropbox/Papers/Murdock Paper/Murdock_Writing/Submissions/Murdock Replication Data")

## SET UP GRID  -----------------------------------------------------------
{
Grid_25deg <- shapefile("Data/Grid Shapefiles/Grid_0_25_degrees.shp")
Grid_5degree <- shapefile("Data/Grid Shapefiles/Grid_5_degree.shp")
Grid_1degree <- shapefile("Data/Grid Shapefiles/Grid_1_degree.shp")
Grid_2degree <- shapefile("Data/Grid Shapefiles/Grid_2_degree.shp")

Grid_1degree_2 <- intersect(Grid_1degree, Grid_25deg)
Grid_1degree_2$Grid_1deg_Area_km <- area(Grid_1degree_2) /1000000
DF <- as.data.frame(Grid_1degree_2)
DF <- DF[ c(6, 8, 9) ]
write.table(DF, file="Appendix/Larger Grids/Grid_1_degree_MERGE_25km.csv", na = "", sep=",",row.names=F)

Grid_2degree_2 <- intersect(Grid_2degree, Grid_25deg)
Grid_2degree_2$Grid_2deg_Area_km <- area(Grid_2degree_2) /1000000
DF <- as.data.frame(Grid_2degree_2)
DF <- DF[ c(6, 8, 9) ]
write.table(DF, file="Appendix/Larger Grids/Grid_2_degree_MERGEE_25km.csv", na = "", sep=",",row.names=F)

Grid_5degree_2 <- intersect(Grid_5degree, Grid_25deg)
Grid_5degree_2$Grid_5deg_Area_km <- area(Grid_5degree_2) /1000000
DF <- as.data.frame(Grid_5degree_2)
DF <- DF[ c(6, 8, 9) ]
write.table(DF, file="Appendix/Larger Grids/Grid_0_5_degree_MERGEE_25km.csv", na = "", sep=",",row.names=F)

Grid_25deg$centroids <- getSpPPolygonsLabptSlots(Grid_25deg)
Grid_25deg$Grid_area_km <- area(Grid_25deg) /1000000
Country <- shapefile("Data/Grid Covariates/Admin Units/SampledCountries.shp")

Centroid_num <- Grid_25deg$centroids
centroid.sp<-Centroid_num
colnames(centroid.sp) <- c("longitude", "latitude")
df <- data.frame(centroid.sp)
centroid_points <- st_as_sf(df, coords = c("longitude","latitude"), crs = 4326)   

Country$Country_Area_km <- area(Country) /1000000
Country2 <- intersect(Country, Grid_25deg)
Country2$Atlas_CountryArea_Intersect <- area(Country2)/1000000

require(data.table) 
group <- as.data.table(Country2)

sorted_group <- group[order(group$Grid_ID, -group$Atlas_CountryArea_Intersect),] 
ans <- sorted_group[!duplicated(sorted_group$Grid_ID),] 

Grid_25deg <- merge(Grid_25deg, ans, by = "Grid_ID")
}

## IMPORT POLITY ESTIMATES -----------------------------------------------------------
{
#Atlas data - all polities
  Atlas_6hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Atlas_6hr.shp")
  Atlas_8hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Atlas_8hr.shp")
  Atlas_10hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Atlas_10hr.shp")
  Atlas_6hr$Atlas_6hr_area_km <- area(Atlas_6hr) /1000000
  names(Atlas_6hr)[names(Atlas_6hr) == "Polity"] <- "Polity_6hr_Clip"
  Atlas_8hr$Atlas_8hr_area_km <- area(Atlas_8hr) /1000000
  names(Atlas_8hr)[names(Atlas_8hr) == "Polity"] <- "Polity_8hr_Clip"
  Atlas_10hr$Atlas_10hr_area_km <- area(Atlas_10hr) /1000000
  names(Atlas_10hr)[names(Atlas_10hr) == "Polity"] <- "Polity_10hr_Clip"
  
#Atlas data -  composite polities (Sokoto, Adamawa, Mossi, Yoruba)
  Atlas_Comp_6hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities/Atlas_Composite_6hr.shp")
  Atlas_Comp_8hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities/Atlas_Composite_8hr.shp")
  Atlas_Comp_10hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities/Atlas_Composite_10hr.shp")
  Atlas_Comp_6hr$Atlas_Comp_6hr_area_km <- area(Atlas_Comp_6hr) /1000000
  names(Atlas_Comp_6hr)[names(Atlas_Comp_6hr) == "Polity"] <- "Polity_Comp_6hr_Clip"
  Atlas_Comp_8hr$Atlas_Comp_8hr_area_km <- area(Atlas_Comp_8hr) /1000000
  names(Atlas_Comp_8hr)[names(Atlas_Comp_8hr) == "Polity"] <- "Polity_Comp_8hr_Clip"
  Atlas_Comp_10hr$Atlas_Comp_10hr_area_km <- area(Atlas_Comp_10hr) /1000000
  names(Atlas_Comp_10hr)[names(Atlas_Comp_10hr) == "Polity"] <- "Polity_Comp_10hr_Clip"
  
#Atlas data - Polities late 1800s
  Atlas_1880_6hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Polities_late1800s/Atlas_1880_6hr.shp")
  Atlas_1880_8hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Polities_late1800s/Atlas_1880_8hr.shp")
  Atlas_1880_10hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Polities_late1800s/Atlas_1880_10hr.shp")
  Atlas_1880_6hr$Atlas_1880_6hr_area_km <- area(Atlas_1880_6hr) /1000000
  names(Atlas_1880_6hr)[names(Atlas_1880_6hr) == "Polity"] <- "Polity_1880_6hr_Clip"
  Atlas_1880_8hr$Atlas_1880_8hr_area_km <- area(Atlas_1880_8hr) /1000000
  names(Atlas_1880_8hr)[names(Atlas_1880_8hr) == "Polity"] <- "Polity_1880_8hr_Clip"
  Atlas_1880_10hr$Atlas_1880_10hr_area_km <- area(Atlas_1880_10hr) /1000000
  names(Atlas_1880_10hr)[names(Atlas_1880_10hr) == "Polity"] <- "Polity_1880_10hr_Clip"
  
#Atlas data -Polities in late 1800s with composite polities (Sokoto, Adamawa, Mossi, Maravi)
  Atlas_Comp_1880_6hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities_late1800s/Atlas_Composite_1880_6hr.shp")
  Atlas_Comp_1880_8hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities_late1800s/Atlas_Composite_1880_8hr.shp")
  Atlas_Comp_1880_10hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities_late1800s/Atlas_Composite_1880_10hr.shp")
  Atlas_Comp_1880_6hr$Atlas_Comp_1880_6hr_area_km <- area(Atlas_Comp_1880_6hr) /1000000
  names(Atlas_Comp_1880_6hr)[names(Atlas_Comp_1880_6hr) == "Polity"] <- "Polity_Comp_1880_6hr_Clip"
  Atlas_Comp_1880_8hr$Atlas_Comp_1880_8hr_area_km <- area(Atlas_Comp_1880_8hr) /1000000
  names(Atlas_Comp_1880_8hr)[names(Atlas_Comp_1880_8hr) == "Polity"] <- "Polity_Comp_1880_8hr_Clip"
  Atlas_Comp_1880_10hr$Atlas_Comp_1880_10hr_area_km <- area(Atlas_Comp_1880_10hr) /1000000
  names(Atlas_Comp_1880_10hr)[names(Atlas_Comp_1880_10hr) == "Polity"] <- "Polity_Comp_1880_10hr_Clip"
}

## IMPORT DVs  -----------------------------------------------------------
Nightlight_2013 <- raster("Data/Outcome Variables/Nightlight/F182013.v4/F182013.v4c_web.avg_vis.tif")
Nightlight_1992 <- raster("Data/Outcome Variables/Nightlight/F101992.v4/F101992.v4b_web.avg_vis.tif")
EWI <- shapefile("Data/Outcome Variables/EstimatedWealth/EstimatedWealthData.shp")


## IMPORT COVARIATES -----------------------------------------------------------
{
  Diamonds <- shapefile("Data/Grid Covariates/Diamonds/DIADATA Data/DIADATA.shp")
  Gems <- shapefile("Data/Grid Covariates/Gem Data/GEMDATA.shp")
  Copper <- shapefile("Data/Grid Covariates/USGS Minerals/USGS_Copper.shp")
  Iron <- shapefile("Data/Grid Covariates/USGS Minerals/USGS_IRon.shp")
  Tin <- shapefile("Data/Grid Covariates/USGS Minerals/USGS_Tin.shp")
  Gold <- shapefile("Data/Grid Covariates/USGS Minerals/USGS_Gold.shp")
  Silver <- shapefile("Data/Grid Covariates/USGS Minerals/USGS_Silver.shp")
  Coastline <- shapefile("Data/Grid Covariates/Water/Africa_Coastline2.shp")
  Water <- shapefile("Data/Grid Covariates/Water/Africa_waterbody.shp")
  MajorRivers <- shapefile("Data/Grid Covariates/Water/mrb_rivers.shp")
  AfricanCapitals <- shapefile("Data/Grid Covariates/Admin Units/African_Capitals.shp")
  AfricanBorders <- shapefile("Data/Grid Covariates/Admin Units/African Borders2.shp")
  SlaveExports_Nunn <- shapefile("Data/Grid Covariates/Nunn/Murdock_Nunn Merge.shp")  
  Water2 <- intersect(Water, Grid_25deg)
  Water2$Water_Area_km <- area(Water2) /1000000
  DF <- as.data.frame(Water2)
  df2 <- DF %>% group_by(Grid_ID) %>% summarize(Water_Area_km2 = sum(Water_Area_km))
  require(sp)
  Grid_25deg <- merge(Grid_25deg,df2, by="Grid_ID")
  Grid_25deg$Grid_PercH20 <- (Grid_25deg$Water_Area_km2/Grid_25deg$Grid_area_km.x)*100
  
  LandSuitability <- nc_open("Data/Grid Covariates/Land Suitability/raw_data_land_suit_land_suit_0.50x0.50.nc")
  print(LandSuitability)
  LandSuit_Lon <- ncvar_get(LandSuitability, "longitude")
  LandSuit_Lat <- ncvar_get(LandSuitability, "latitude", verbose = F)
  LandSuit_t <- ncvar_get(LandSuitability, "time")
  LandSuit <- ncvar_get(LandSuitability, "data")
  fillvalue <- ncatt_get(LandSuitability, "data", "_FillValue")
  fillvalue
  LandSuit[LandSuit == fillvalue$value] <- NA
  LandSuitability <- raster(t(LandSuit), xmn=min(LandSuit_Lon), xmx=max(LandSuit_Lon), ymn=min(LandSuit_Lat), ymx=max(LandSuit_Lat), crs = 4326)
  LandSuitability_2 <- crop(LandSuitability, Grid_25deg, snap="out")

  Precipitation <- nc_open("Data/Grid Covariates/Precipitation/chirps-v2.0.annual.nc")
  Precipitation_lon <- ncvar_get(Precipitation, "longitude")
  Precipitation_lat <- ncvar_get(Precipitation, "latitude", verbose = F)
  Prec.array <- ncvar_get(Precipitation, "precip") 
  dim(Prec.array) 
  fillvalue <- ncatt_get(Precipitation, "precip", "_FillValue")
  fillvalue
  Prec.array[Prec.array == fillvalue$value] <- NA
  Prec.array <- Prec.array[, , 1] 
  Precipitation <- raster(t(Prec.array), xmn=min(Precipitation_lon), xmx=max(Precipitation_lon), ymn=min(Precipitation_lat), ymx=max(Precipitation_lat), crs = 4326)
  Precipitation2  <- terra::flip(Precipitation, direction="y")
  Precipitation_2 <- crop(Precipitation2, Grid_25deg, snap="out")
}

## SAMPLE BY GRID  -----------------------------------------------------------
# Polity estimates
{
  Atlas_6hr_Intersect <- intersect(Atlas_6hr, Grid_25deg)
  Atlas_8hr_Intersect <- intersect(Atlas_8hr, Grid_25deg)
  Atlas_10hr_Intersect <- intersect(Atlas_10hr, Grid_25deg)
  Atlas_6hr_Intersect$Atlas_6hr_area_km_Intersect <- area(Atlas_6hr_Intersect) /1000000
  Atlas_8hr_Intersect$Atlas_8hr_area_km_Intersect <- area(Atlas_8hr_Intersect) /1000000
  Atlas_10hr_Intersect$Atlas_10hr_area_km_Intersect <- area(Atlas_10hr_Intersect) /1000000
  
  Atlas_Comp_6hr_Intersect <- intersect(Atlas_Comp_6hr, Grid_25deg)
  Atlas_Comp_8hr_Intersect <- intersect(Atlas_Comp_8hr, Grid_25deg)
  Atlas_Comp_10hr_Intersect <- intersect(Atlas_Comp_10hr, Grid_25deg)
  Atlas_Comp_6hr_Intersect$Atlas_Comp_6hr_area_km_Intersect <- area(Atlas_Comp_6hr_Intersect) /1000000
  Atlas_Comp_8hr_Intersect$Atlas_Comp_8hr_area_km_Intersect <- area(Atlas_Comp_8hr_Intersect) /1000000
  Atlas_Comp_10hr_Intersect$Atlas_Comp_10hr_area_km_Intersect <- area(Atlas_Comp_10hr_Intersect) /1000000
  
  Atlas_1880_6hr_Intersect <- intersect(Atlas_1880_6hr, Grid_25deg)
  Atlas_1880_8hr_Intersect <- intersect(Atlas_1880_8hr, Grid_25deg)
  Atlas_1880_10hr_Intersect<- intersect(Atlas_1880_10hr, Grid_25deg)
  Atlas_1880_6hr_Intersect$Atlas_1880_6hr_area_km_Intersect <- area(Atlas_1880_6hr_Intersect) /1000000
  Atlas_1880_8hr_Intersect$Atlas_1880_8hr_area_km_Intersect <- area(Atlas_1880_8hr_Intersect) /1000000
  Atlas_1880_10hr_Intersect$Atlas_1880_10hr_area_km_Intersect <- area(Atlas_1880_10hr_Intersect) /1000000
  
  Atlas_Comp_1880_6hr_Intersect <- intersect(Atlas_Comp_1880_6hr, Grid_25deg)
  Atlas_Comp_1880_8hr_Intersect <- intersect(Atlas_Comp_1880_8hr, Grid_25deg)
  Atlas_Comp_1880_10hr_Intersect <- intersect(Atlas_Comp_1880_10hr, Grid_25deg)
  Atlas_Comp_1880_6hr_Intersect$Atlas_Comp_1880_6hr_area_km_Intersect <- area(Atlas_Comp_1880_6hr_Intersect) /1000000
  Atlas_Comp_1880_8hr_Intersect$Atlas_Comp_1880_8hr_area_km_Intersect <- area(Atlas_Comp_1880_8hr_Intersect) /1000000
  Atlas_Comp_1880_10hr_Intersect$Atlas_Comp_1880_10hr_area_km_Intersect <- area(Atlas_Comp_1880_10hr_Intersect) /1000000
}

# Outcome Variables
  Grid_25deg$Nightlight_1992_mean <- exact_extract(Nightlight_1992,Grid_25deg,fun="mean")
  Grid_25deg$Nightlight_20132_mean <- exact_extract(Nightlight_2013,Grid_25deg,fun="mean")
  EWI_Intersect <- intersect(EWI, Grid_25deg)

# Covariates
{
  Grid_25deg$Diamonds <- over(SpatialPolygons(Grid_25deg@polygons), SpatialPoints(Diamonds))
  Grid_25deg$Gems <- over(SpatialPolygons(Grid_25deg@polygons), SpatialPoints(Gems))
  Grid_25deg$Silver <- over(SpatialPolygons(Grid_25deg@polygons), SpatialPoints(Silver))
  Grid_25deg$Copper <- over(SpatialPolygons(Grid_25deg@polygons), SpatialPoints(Copper))
  Grid_25deg$Tin <- over(SpatialPolygons(Grid_25deg@polygons), SpatialPoints(Tin))
  Grid_25deg$Iron <- over(SpatialPolygons(Grid_25deg@polygons), SpatialPoints(Iron))
  Grid_25deg$Gold <- over(SpatialPolygons(Grid_25deg@polygons), SpatialPoints(Gold))

  Mountains_Binary <- raster("Data/Grid Covariates/Global Mountain Explorer/GlobalMountainsK3Binary/k3binary.tif")
  Mountains_Binary2 <- crop(Mountains_Binary, Grid_25deg, snap="out")
  Mountains_Type <- raster("Data/Grid Covariates/Global Mountain Explorer/GlobalMountainsK3Classes/k3classes.tif")
  Mountains_Type2 <- crop(Mountains_Type, Grid_25deg, snap="out")
  Malaria_2000 <- raster("Data/Grid Covariates/Malaria Area Project/202206_Global_Pf_Parasite_Rate_2000/202206_Global_Pf_Parasite_Rate_2000.tif")
  Malaria_2000_2 <- crop(Malaria_2000, Grid_25deg, snap="out")

  PopDensity_0AD <- raster("Data/Grid Covariates/HYDE/Hyde3_2/0AD_pop/popd_0AD.asc")
  PopDensity_0AD_2 <- crop(PopDensity_0AD, Grid_25deg, snap="out")
  
  Elevation <- raster("Data/Grid Covariates/Elevation/Africa_dem_15s.tif")
  Ruggedness <- raster("Data/Grid Covariates/Ruggedness/tri.tif")
  Ruggedness_MeanSlope <- raster("Data/Grid Covariates/Ruggedness/slope.tif")
  Ruggedness_CellArea <- raster("Data/Grid Covariates/Ruggedness/cellarea.tif")
  
  Grid_25deg$mt_binary_mean <- exact_extract(Mountains_Binary2,Grid_25deg,fun="mean")
  Grid_25deg$Grid_MtType <- exact_extract(Mountains_Type2,Grid_25deg,fun="mean")
  Grid_25deg$Malaria_2000 <- exact_extract(Malaria_2000_2,Grid_25deg,fun="mean")
  Grid_25deg$PopDensity_0AD <- exact_extract(PopDensity_0AD_2,Grid_25deg,fun="mean")
  Grid_25deg$Elevation <- exact_extract(Elevation,Grid_25deg,fun="mean")
  Grid_25deg$LandSuitability <- exact_extract(LandSuitability_2,Grid_25deg,fun="mean")
  Grid_25deg$Precipitation <- exact_extract(Precipitation_2,Grid_25deg,fun="mean")
  
  Ruggedness_SlopeWeight<- overlay(Ruggedness_MeanSlope,
                                   Ruggedness_CellArea,
                                   fun=function(r1, r2){return(r1*r2)})
  
  Ruggedness_RuggedWeight<- overlay(Ruggedness,
                                    Ruggedness_CellArea,
                                    fun=function(r1, r2){return(r1*r2)})
  
  Grid_25deg$Rugged_CellAreaTtl <- exact_extract(Ruggedness_CellArea,Grid_25deg,fun="sum")
  
  Grid_25deg$Rugged_SlopeTtl <- exact_extract(Ruggedness_SlopeWeight,Grid_25deg,fun="sum")
  Grid_25deg$Rugged_TriTtl <- exact_extract(Ruggedness_RuggedWeight,Grid_25deg,fun="sum")
  Grid_25deg$CellArea <- exact_extract(Ruggedness_CellArea,Grid_25deg,fun="sum")
  
  ## NEXT CALC GRID AREA + WATER AREA
  Ruggedness2 <- crop(Ruggedness, Grid_25deg, snap="out")
  Ruggedness_CellArea2 <- crop(Ruggedness_CellArea, Grid_25deg, snap="out")
  Ruggedness_MeanSlope2 <- crop(Ruggedness_MeanSlope, Grid_25deg, snap="out")
}
  
# Merge & Export  
  Grid_25deg2 = st_as_sf(Grid_25deg)
  SlaveExports_Nunn2 <- st_as_sf(SlaveExports_Nunn)
  EWI_Intersect2 <- st_as_sf(EWI_Intersect)
  
  centroids2 = st_as_sf(centroid_points)
  Coastline2 = st_as_sf(Coastline)
  centroids2$d_coast_m <- st_distance(x = centroids2, y = Coastline2)
  
  Grid_25deg2 <- st_join(Grid_25deg2, centroids2, join = st_intersects)
  Grid_25deg2 <- st_join(Grid_25deg2, SlaveExports_Nunn2, join = st_intersects)
  Grid_25deg2 <- st_join(Grid_25deg2, EWI_Intersect2, join = st_intersects)
  
  sf_use_s2(FALSE)
  
  Grid_25deg2 <- st_join(Grid_25deg2, centroid_points, join = st_contains)
  
save(Grid_25deg2, file = "Appendix/Larger Grids/GridSampling_PrePolity_25degrees.RData") 


### TAKE THIS DATA AND MERGE IN SETS OF POLITIES BECAUSE WILL PROLIFERATE IF OVERLAP - ESP IN BASE
# Atlas
{
  load("Appendix/Larger Grids/GridSampling_PrePolity_25degrees.RData")

  sf_use_s2(FALSE)
  
  Atlas_8hr_Intersect2 = st_as_sf(Atlas_8hr_Intersect)

  Atlas_8hr_Intersect2 <- Atlas_8hr_Intersect2 %>% select(Polity_8hr_Clip, Atlas_8hr_area_km_Intersect)

  Grid_25deg2 <- st_join(Grid_25deg2, Atlas_8hr_Intersect2, join = st_intersects)

  d <- as.data.frame(Grid_25deg2) 
  d <- d[ -c(75) ]
  write.csv(d, "Appendix/Larger Grids/Atlas_25degrees.csv", row.names = FALSE, na = "")
}

### IMPORT VIOLENCE DATA & MERGE TO GRID
{
  #ACLED data
  ACLED_Base <- 
    read.csv("Data/Outcome Variables/Conflict Data/1997-01-01-2023-03-01-East_Asia-Eastern_Africa-Middle_Africa-Southern_Africa-Western_Africa_MW.csv",
             stringsAsFactors = FALSE)
  
  ACLED_Base_Plotted <- SpatialPointsDataFrame(ACLED_Base[,1:2], ACLED_Base)
  crs(ACLED_Base_Plotted) <- CRS('+init=EPSG:4326')                                   
  #plot(ACLED_Base_Plotted)
  
  ACLED_Base_Plotted2 = st_as_sf(ACLED_Base_Plotted)
  
  ACLED_Base_Plotted2 %>% select(33) -> ACLED_Base_Plotted2
  
  Grid_25deg2_1 <- st_join(Grid_25deg2, ACLED_Base_Plotted2, join = st_intersects)
  
  d <- as.data.frame(Grid_25deg2_1) 
  d <- d[ -c(73) ]
  write.csv(d, "Appendix/Larger Grids/AcledID_25degrees.csv", row.names = FALSE, na = "")
  
  
  #UCDP GED data
  GED_Base <- 
    read.csv("Data/Outcome Variables/Conflict Data/GED_23.1_Africa_Final_MW.csv",
             stringsAsFactors = FALSE)
  GED_Base_Plotted <- SpatialPointsDataFrame(GED_Base[,1:2], GED_Base)
  crs(GED_Base_Plotted) <- CRS('+init=EPSG:4326')                                   
  #plot(GED_Base_Plotted)
  
  GED_Base_Plotted2 = st_as_sf(GED_Base_Plotted)
  
  GED_Base_Plotted2 %>% select(3) -> GED_Base_Plotted2
  
   year <- 1966
  regcap1966 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                       band =  year-1945+1)
  natcap1966 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                       band =  year-1945+1)
  year <- 1977
  regcap1977 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                       band =  year-1945+1)
  natcap1977 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                       band =  year-1945+1)
  
  year <- 1980
  regcap1980 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                       band =  year-1945+1)
  natcap1980 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                       band =  year-1945+1)
  year <- 1968
  regcap1968 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                       band =  year-1945+1)
  natcap1968 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                       band =  year-1945+1)
  year <- 1975
  regcap1975 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                       band =  year-1945+1)
  natcap1975 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                       band =  year-1945+1)
  
  Grid_25deg2$regcap1975 <- exact_extract(regcap1975,Grid_25deg2,fun="mean")
  Grid_25deg2$natcap1975 <- exact_extract(natcap1975,Grid_25deg2,fun="mean")
  Grid_25deg2$regcap1968 <- exact_extract(regcap1968,Grid_25deg2,fun="mean")
  Grid_25deg2$natcap1968 <- exact_extract(natcap1968,Grid_25deg2,fun="mean")
  Grid_25deg2$regcap1980 <- exact_extract(regcap1980,Grid_25deg2,fun="mean")
  Grid_25deg2$natcap1980 <- exact_extract(natcap1980,Grid_25deg2,fun="mean")
  Grid_25deg2$regcap1977 <- exact_extract(regcap1977,Grid_25deg2,fun="mean")
  Grid_25deg2$natcap1977 <- exact_extract(natcap1977,Grid_25deg2,fun="mean")
  Grid_25deg2$regcap1966 <- exact_extract(regcap1966,Grid_25deg2,fun="mean")
  Grid_25deg2$natcap1966 <- exact_extract(natcap1966,Grid_25deg2,fun="mean")
  
  Grid_25deg2_2 <- st_join(Grid_25deg2, GED_Base_Plotted2, join = st_intersects)
  
  sf_use_s2(FALSE)
  Grid_25deg2_2 = st_as_sf(Grid_25deg2)
  
  d <- as.data.frame(Grid_25deg2_2) 
  d <- d[ -c(73) ]
  write.csv(d, "Appendix/Larger Grids/GEDID_25degrees.csv", row.names = FALSE, na = "")
}

#META data
{
  Meta_Base <- 
    read.csv("Data/Outcome Variables/Meta Relative Wealth Data/Africa_Master_relative_wealth_index.csv",
             stringsAsFactors = FALSE)
  Meta_Base_Plotted <- SpatialPointsDataFrame(Meta_Base[,3:2], Meta_Base)
  crs(Meta_Base_Plotted) <- CRS('+init=EPSG:4326')                                   
  plot(Meta_Base_Plotted)
  
  sf_use_s2(FALSE)
  
  Meta_Base_Plotted2 = st_as_sf(Meta_Base_Plotted)
  
  Grid_25deg2_2 <- st_join(Grid_25deg2, Meta_Base_Plotted2, join = st_intersects)
  
  d <- as.data.frame(Grid_25deg2_2) 
  d <- d[ -c(36) ]
  write.csv(d, "Appendix/Larger Grids/Meta_Wealth_25degrees.csv", row.names = FALSE, na = "")
}



